#####################################################################
# Calibration and DGPs based on Andersson (2019)
#####################################################################

# Packages

library(foreign)
library(pracma)

# Data

carbontax_data <- read.dta("data/carbontax_data.dta")

keeps <- c("Countryno", "year","CO2_transport_capita")
carbontax_data <- carbontax_data[keeps]

carbontax_wide <- reshape(carbontax_data, timevar = "year",idvar = c("Countryno"),direction = "wide")

Y1_app <- as.matrix(cbind(carbontax_wide[13,-1]))
Y0_app <- t(as.matrix(carbontax_wide[-13,-1]))

T0_app  <- 30 # 1960-1989
T1_app  <- length(Y1_app)-T0_app
T01_app <- T0_app + T1_app

J_app <- dim(Y0_app)[2]

# Calibration

w0_sc         <- as.matrix(scinference:::sc(Y1_app[1:T0_app],Y0_app[1:T0_app,],lsei_type = 1)$w.hat)
u_hat         <- Y1_app[1:T0_app] - Y0_app[1:T0_app,] %*% w0_sc
ar_obj_u_hat  <- ar(u_hat, order.max=1)
rho_u         <- ar_obj_u_hat$ar 
var_u         <- ar_obj_u_hat$var.pred 

Y0dt  <- apply(Y0_app,2,detrend)

result <- svd(Y0dt)
max(abs(result$u %*% diag(result$d) %*% t(result$v)-Y0dt))

num_factor            <- 4
dd                    <- result$d
dd[(num_factor+1):length(result$d)]  <- 0
Y0_PCA_est            <- result$u %*% diag(dd) %*% t(result$v)
Y0_resid              <- Y0dt-Y0_PCA_est

Factors <- result$u[,1:num_factor]
Lambda  <- result$v %*% diag(dd)
Lambda  <- Lambda[,1:num_factor]

var_factors <- matrix(NA,num_factor,1)
for (j in 1:num_factor){
  var_factors[j]  <- var(Factors[,j])
}

rho_vec       <- matrix(0,J_app,1) # default is =0, and it gets replaced only if we find a positive AR(1) coefficient
var_epsl_vec  <- matrix(NA,J_app,1)
for (j in 1:J_app){
  ar_obj  <- ar(Y0_resid[,j], order.max=1)
  if (ar_obj$order>0){rho_vec[j,1] <- ar_obj$ar}
  var_epsl_vec[j] <- ar_obj$var.pred
}

# Print output (commented out to save space)

# print("SC weights")
# print(round(w0_sc,digits=2))
# print("AR(1) coefficient treated unit")
# print(round(rho_u,digits=2))
# print("Summary statistics of AR(1) coefficients for control units")
# print(round(c(min(rho_vec),max(rho_vec),median(rho_vec)),digits=2))
# 


